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Abstract. We present results of the all-sky search for gravitational-wave signals 
from spinning neutron stars in the data of the EXPLORER resonant bar detector. 
Our data analysis technique was based on the maximum likelihood detection 
method. We briefly describe the theoretical methods that we used in our 
search. The main result of our analysis is an upper limit of 2 X 10~^^ for the 
dimensionless amplitude of the continuous gravitational- wave signals coming from 
any direction in the sky and in the narrow frequency band from 921.00 Hz to 
921.76 Hz. 
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1. Introduction 

A unique property of the gravitational-wave detectors is that with a single observation 
resulting in one time series one can search for signals coming from all sky directions. 
In the case of other instruments like optical and radio telescopes, to cover the whole 
or even part of the sky requires a large amount of expensive telescope time as each 
sky position needs to be observed independently. The difficulty in the search for 
gravitational- wave signals is that they are very weak and they are deeply buried in the 
noise of the detector. Consequently the detection of these signals and interpretation 
of data analysis results is a delicate task. In this paper we present results of an 
all-sky search for continuous sources of gravitational radiation. A prime example 
of such a source is a spinning neutron star. A signal form such a source has definite 
characteristics that make it suitable for application of the optimal detection techniques 
based on matched filtering. Moreover such signals are stable as a result of the stability 
of the rotation of the neutron star and they will be present in the data for time periods 
much longer than the observational interval. This enables a reliable verification of the 
potential candidates by repeating the observations both by the same detector and 
by different detectors. To perform our all-sky search we have used the data of the 
EXPLORER J] resonant bar detector. The directed search of the galactic center 
with the EXPLORER detector has already been carried out and an upper limit for 
the amplitude of the gravitational waves has been established 2, . 

Our paper is divided into two parts. In the first part we summarize the theoretical 
tools that we use in our analysis and in the second part we present results of our all-sky 
search. The main result of our analysis is an upper limit of 2 x 10~^^ for dimensionless 
amplitude of gravitational waves originating from continuous sources located in any 
position in the sky and in the frequency band from 921.00 Hz to 921.76 Hz. 

The data analysis was performed by a team consisting of Pia Astone, 
Kaz Borkowski, Piotr Jaranowski, Andrzej Krolak and Maciej Pietka and 
was carried out on the basis of Memorandum of Understanding between 
the ROG collaboration and Institute of Mathematics of Polish Academy of 
Sciences. More details about the search can be found on the website: 
|http : //www, astro .uni . toruriT pl/~kb/AllSky/AllSky . htinl| 

2. Data analysis methods 

In this section we give a summary of data analysis techniques that we used in the 
search. The full exposition of our method is given in Ref. 

2.1. Response of a bar detector to a continuous gravitational-wave signal 

Dimensionless noise-free response function /i of a resonant bar gravitational-wave 
detector to a weak plane gravitational wave in the long wavelength approximation 
[i.e., when the size of the detector is much smaller than the reduced wavelength A/ (27r) 
of the wave] can be written as a linear combination of the wave polarization functions 
h+ and hx ■ 



where and -Fx are called the beam-pattern functions and can be written as 



h{t)^F+{t) h+{t)+Fx{t)hx{t), 



(1) 



a{t) cos 2x1; -f b{t) sin 2?/;, 
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Fy,{t) = b{t) cos 2ij - a{t) sin 2^Ij, (2) 
where ip is the polarization angle of the wave and 

a{t) = ^(cos^ 7 - sin^ 7sin^ (/}){! + sin^ S) cos[2(a ~ 4ir - flrt)] 
+ ^ sin 27 sin 0(1 + sin'^ S) sin[2(a — (j)r — Qrt)] 

— ^ sin^ 7 sin 20 sin 2(5 cos(q: — 0^ — f^r^) 
+ - sin 27 cos sin 28 sin (a — </>,. — fi^i) 

+ - ( 1 — 3 sin^ 7 cos^ 0) cos^ (5, 

h(t) = — sin27sin(/)sin(5cos[2(a — 0^ — ^^rO] 

+ (cos^ 7 — sin^ 7sin^ 0) sin(5sin[2(a — 4>r — ^A)] 

— sin 27 cos cos 5 cos(a — <j)r — flrt) 

— sin^ 7 sin 20 cos (5 sin(a — 0^ — fij,^) ■ (3) 

In Eqs. lO and (|3Jl the angles a and ^ are respectively right ascension and declination 
of the gravitational-wave source. The geodetic latitude of the detector's site is denoted 
by 0, the angle 7 determines the orientation of the bar detector with respect to local 
geographical directions: 7 is measured counter-clockwise from East to the bar's axis 
of symmetry. The rotational angular velocity of the Earth is denoted by Q^, f^nd 0,, 
is a deterministic phase which defines the position of the Earth in its diurnal motion 
at t = (the sum 0^ + ilrt essentially coincides with the local sidereal time of the 
detector's site, i.e., with the angle between the local meridian and the vernal point). 

We are interested in a continuous gravitational wave, which is described by the 
wave polarization functions of the form 

h+{t) = ha+cos^(t), /ix(<) = /lox sin^'(t), (4) 

where /10+ and /iqx are independent constant amplitudes of the two wave polarizations. 
These amplitudes depend on the physical mechanisms generating the gravitational 
wave. In the case of a wave originating from a spinning neutron star these amplitudes 
can be estimated by 




where / is the neutron star moment of inertia with respect to the rotation axis, e is 
the poloidal ellipticity of the star, and r is the distance to the star. The value of 10~^ 
of the parameter e in the above estimate should be treated as an upper bound. In 
reality it may be several orders of magnitude less. 
The phase 5' of the wave is given by 

*(t) = $o + $(i), 
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joining the solar system barycenter (SSB) with the detector, no is the constant unit 
vector in the direction from the SSB to the source of the gravitational-wave signal. 
We assume that the gravitational-wave signal is almost monochromatic around some 
angular frequency ujq which we define as instantaneous angular frequency evaluated 
at the SSB a,t t — 0, oj^ (fc = 1, 2, . . .) is the fcth time derivative of the instantaneous 
angular frequency at the SSB evaluated at t = 0. To obtain formulas © we model 
the frequency of the signal in the rest frame of the neutron star by a Taylor series. 
For the detailed derivation of the phase model © see Sec. II B and Appendix A of 
Ref. [H]. The integers si and S2 are the number of spin downs to be included in the 
two components of the phase. We need to include enough spin downs so that we have 
a sufficiently accurate model of the signal to extract it from the noise. 

The vector rssB is the sum of two vectors: the vector tes joining the solar system 
barycenter and the Earth barycenter and vector te joining the Earth barycenter and 
the detector. The vector Tes can be obtained form the JPL Planetary and Lunar 
Ephemerides: "DE405/LE405" . They are available via the Internet at anonymous 
ftp: navigator, jpl.nasa.gov, the directory: ephem/export. Whereas the vector Te 
can be accurately computed using the International Earth Rotation Service (lERS) 
tables available at http://hpiers.obspm.fr/iers/eop/eopc04. The codes to read 
the above files and to calculate the vector rsse are described in Ref. [Hj. 

It is convenient to write the response of the gravitational-wave detector given 
above in the following form 



where Ai are four constant amplitudes, the functions a and b are given by Eqs. 
and $ is the phase given by Eq. ©. The modulation amplitudes a and b depend on 
the right ascension a and the declination 5 of the source (they also depend on the 
angles </> and 7). The phase $ depends on the frequency luq, s spin-down parameters 
ujk (fc = 1, . . . , s), and on the angles a, S. We call the parameters too, uJk, a, 6 the 
intrinsic parameters and the remaining ones the extrinsic parameters. Moreover the 
phase $ depends on the latitude cj) of the detector. The whole signal h depends on 
s + 5 unknown parameters: /10+, /iqx, a, S, too, uJk- 

2.2. Optimal data analysis method 

We assume that the noise n in the detector is an additive, stationary, Gaussian, and 
zero-mean random process. Then the data x (when the signal h is present) can be 
written as 



h{t) = Aia{t) cos ^{t) + A2b[t) cos ^(t) + A'ia[t) sin $(i) + A^bit) sin $(i). 



(7) 



x{t) = n{t) + h{t). 
The log likelihood function has the form 



(8) 



\ogK^{x\h)--{h\h), 



(9) 



where the scalar product ( • | ■ ) is defined by 




(10) 



All-sky upper limit 



5 



In A 



x{t)dt. (12) 



In Eq. H1U|) denotes the real part, tilde denotes Fourier transform, asterisk is complex 
conjugation, and Sh is the one-sided spectral density of the detector's noise. Assuming 
that over the bandwidth of the signal h the spectral density Sh{f) is nearly constant 
and equal to Sh{fo)i where /o is the frequency of the signal measured at the SSB at 
t = 0, the log likelihood ratio from Eq. (j^ can be approximated by 

where (•) denotes time averaging over the observational interval [0,To]: 
1 

J^o Jo 

The signal h depends linearly on four amplitudes Ai. The equations for the 
maximum likelihood (ML) estimators Ai of the amplitudes Ai are given by 

.^1,...,4. ,13, 

One can easily find the explicit analytic solution to Eqs. To simplify formulas 

we assume that the observation time To is an integer multiple of one sidereal day, i.e., 
To = n{2n/flr) for some positive integer n. Then the time average of the product 
of the functions a and b vanishes, (ab) = 0, and the analytic formulas for the ML 
estimators of the amplitudes are given by 

A-W' ^'-^W' ^'-^w- 

The reduced log likelihood function J-' is the log likelihood function where amplitude 
parameters Ai were replaced by their estimators Ai. By virtue of Eqs. (|14|l from Eq. 
l(TT|l one gets 

where 

Fa-.^ [ x{t)a{t)exp[-i<^{t)]dt, 



Fb := [ x{t)b{t)exp[~i^{t)]dt. (16) 
Jo 

The ML estimators of the signal parameters are obtained in two steps. Firstly, 
the estimators of the frequency, the spin-down parameters, and the angles a and 6 are 
obtained by maximizing the functional with respect to these parameters. Secondly, 
the estimators of the amplitudes Ai are calculated from the analytic formulas (|14|l 
with the correlations (xhi) evaluated for the values of the parameters obtained in the 
first step. 

In order to calculate the optimum statistics J- efficiently we introduce an 
approximation to the phase of the signal that is valid for observation times short 
comparing to the period of 1 year. The phase of the gravitational-wave signal contains 
terms arising from the motion of the detector with respect to the SSB. These terms 
consist of two contributions, one which comes from the motion of the Earth barycenter 
with respect to the SSB, and the other which is due to the diurnal motion of the 
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detector with respect to the Earth barycenter. The first contribution has a period of 
one year and thus for shorter observation times can be well approximated by a few 
terms of the Taylor expansion. The second term has a period of 1 sidereal day and 
to a very good accuracy can be approximated by a circular motion. We find that for 
two days of observation time that we used in the analysis of the EXPLORER data 
presented in Sectional we can truncate the expansion at terms that are quadratic in 
time. Moreover in the Taylor expansion of the frequency parameter it is sufficient to 
include terms up to the first spin down. We thus propose the following approximate 
simple model of the phase of the gravitational-wave signal: 

'l's(i) = p + Pot + pit^ + Acos{Vlrt) + Bsm{Q.rt), (17) 

where Vt^. is the rotational angular velocity of the Earth. The characteristic feature 
of the above approximation is that the phase is a linear function of the parameters of 
the signal. The parameters A and B can be related to the right ascension a and the 
declination 5 of the gravitational-wave source through the equations 

A ~ cos 5 cos( a — 4>r) ^ (18) 
c 

i? = cos(5sin(a — (/)r), (19) 
c 

where ujq is the angular frequency of the gravitational-wave signal and r is the 
equatorial component of the detector's radius vector. The parameters p, Pq, and 
Pi contain contributions both form the intrinsic evolution of the gravitational-wave 
source and the modulation of the signal due to the motion of the Earth around the 
Sun. With this approximation the integrals given by Eqs. l(TB|) and needed to compute 
^ become Fourier transforms and they can be efficiently calculated using the EFT 
algorithm. Thus the evaluation of J-' consists of correlation of the data with two linear 
filters depending on parameters pi,A,B followed by FFTs. In Ref. we have verified 
that for the case of our search the linear approximation to the phase does not lead to 
the loss of signal-to-noise ratio of more than 5%. 

In order to identify potential gravitational-wave signals we apply a two step 
procedure consisting of a coarse search followed by a fine search. The coarse search 
consists of evaluation of on a discrete grid constructed in such a way that the 
loss of the signal-to-noise is minimized and comparison of the obtained values of !F 
with a predefined threshold To- The parameters of the nodes of the grid for which 
the threshold is crossed are registered as potential gravitational- wave signals. The 
threshold is calculated from a chosen false alarm probability which is defined as the 
probability that crosses the threshold when no signal is present and the data is only 
noise. The fine search consists of finding a local maximum of using a numerical 
implementation of the Nelder-Mead algorithm, where coordinates of the starting point 
of the maximization procedure are the parameter values of the coarse search. 

To calculate the false alarm probability as a function of a threshold and to 
construct a grid in the parameter space we introduce yet another approximation of 
our signal. Namely we use a signal with a constant amplitude and the phase given 
by Eq. (|17|) . In paper 0] we have shown that the Fisher matrix for the exact model 
with amplitude modulations given by Eq. ^ can be accurately reproduced by the 
Fisher matrix of a constant amplitude model. As the calculations of the false alarm 
probability and construction of a grid in the parameter space are based on the Fisher 
matrix we expect that the constant amplitude model is a good approximation for the 
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purpose of the above calculations. We stress that in the search of real data we used 
the full model with amplitude modulations. 

We shall now summarize the statistical properties of the functional T. We first 
assume that the phase parameters are known and consequently that T is only a 
function of the random data x. When the signal is absent 2T has a distribution 
with four degrees of freedom and when the signal is present it has a noncentral 
distribution with four degrees of freedom and noncentrality parameter equal to the 
optimal signal-to-noise ratio d defined as 

d := ^/{h\h). (20) 

Consequently the probability density functions po and pi when respectively the signal 
is absent and present arc given by 

Po{T) =Tcxp{-T), (21) 

(22) 



Pi{d,T) = l^I^(dVW) exp (^-T - ^d' 



where /i is the modified Bessel function of the first kind and order 1. The false alarm 
probability Pp is the probability that T exceeds a certain threshold To when there is 
no signal. In our case we have 

POO 

Pf{To) := J PoiT) rfJT = (1 + To) exp{-To). (23) 

The probability of detection Pn is the probability that T exceeds the threshold To 
when the signal-to-noise ratio is equal to d: 

POO 

PD{d,To):= pi{d,T)dT. (24) 

The integral in the above formula cannot be evaluated in terms of known special 
functions. 

We see that when the noise in the detector is Gaussian and the phase parameters 
are known the probability of detection of the signal depends on a single quantity: the 
optimal signal-to-noise ratio d. When the phase parameters are unknown the optimal 
statistics T depends not only on the random data x but also on the phase parameters 
that we shall denote by Such an object is called in the theory of stochastic processes 
a random field. Let us consider the correlation function of the random field 

C{^,0 ■.= E{[T{0 - mmm') - miO]} , (25) 

where 

m{0:=E{TiO}. (26) 
For the case of the constant amplitude, linear phase model we have 

m{0 = 1. (27) 

and 

= {cos[Apot + Apit"^ + AAcos{nrt) + AS sin(f7rt)])^ 
+ (sin[Apo t + Api t^ + AA cosiflrt) + AB sin(fir0])t?8) 
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where A denotes the difference between the parameter values. Thus the correlation 
function C depends only on the difference £,—£,' and not on the values of the parameters 
themselves. In this case the random field J- is called homogeneous. For such fields 
we have developed ^ a method to calculate the false alarm probability that we only 
summarize here. The main idea is to divide the space of the intrinsic parameters ^ into 
elementary cells. The size of the cell is determined by the characteristic correlation 
hypersurface of the random field J-. The correlation hypersurface is defined by the 
condition that at this hypersurface the correlation C equals to half its maximum 
value. Assuming that C attains its maximum value when — ^' = the equation of 
the characteristic correlation hypersurface reads 

C{t) = ^C{0), (29) 

where we have introduced r :— ^ — . Let us expand the left-hand side of Eq. (|29|l 
around t = up to terms of the second order in r. We arrive at the equation 

M 

G^JnT, = 1, (30) 

where M is the dimension of the intrinsic parameter space and the matrix G is defined 
as follows 

1 d'Cir) 
■ C(0) dndr, 



(31) 

T = 

Equation (|30|l defines the boundary of an M-dimensional hyperellipsoid which we 
call the correlation hyperellipsoid. The M-dimensional hypervolume Vceii of the 
hyperellipsoid defined by Eq. (|3U|) equals 

Vcell = ; , (32) 

r(Af/2 + l)VditG ^ ' 

where F denotes the Gamma function. We estimate the number Nc of elementary cells 
by dividing the total Euclidean volume Vtotai of the parameter space by the volume 
Vcell of the correlation hyperellipsoid, i.e. we have 

Nc = -5^. (33) 

Vcell 

We approximate the probability distribution of J-{^) in each cell by the probability 
distribution po{^) when the parameters are known [in our case it is the one given by 
Eq. (|^ ] . The values of the statistics in each cell can be considered as independent 
random variables. The probability that does not exceed the threshold To in a given 
cell is 1 — Pf{To)i where Pf{J^o) is given by Eq. Consequently the probability 

that T does not exceed the threshold To in all the Nc cells is [1 — PpiTo)]^" . The 
probability Pp that T exceeds To in one or more cells is thus given by 

P^{To) = l-[l-PF{To)f^. (34) 

This is the false alarm probability when the phase parameters are unknown. 

The basic quantity to consider in the construction of the grid of templates in the 
parameter space is the expectation value Ei {T} of the statistics when the signal is 
present. For the constant amplitude linear phase model we have that 

El {T{0} - 1 + ^d'C(0, (35) 
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where C is the correlation function of the random field introduced above and d is the 
signal-to-noise ratio of the signal present in the data. We choose the grid of templates 
in such a way that the correlation between any potential signal present in the data 
and the nearest point of the grid never falls below a certain value. In the case of the 
approximate model of the signal that we use the grid is unform and consists of regular 
polygons in the space parametrized by po,pi. A, B. The construction of the grid is 
described in detail in Section VIIA of 0. For the grid that we used in the search of 
the EXPLORER data the correlation function C for any signal present in the data 
was greater than 0.77. 

3. An all-sky search of the EXPLORER data 

We have implemented the theoretical tools presented in Section |21 and we have 
performed an all-sky search for continuous sources of gravitational waves in the data of 
the resonant bar detector EXPLORERf [T]. The EXPLORER detector has collected 
many years of data with a high duty cycle (e.g. in 1991 the duty cycle was 75%). The 
EXPLORER detector was most sensitive over certain two narrow bandwidths (called 
minus and plus modes) of about 1 Hz wide at two frequencies around 1 kHz. To make 
the search computationally manageable we analyzed two days of data in the narrow 
band were the detector had the best sensitivity. By narrowing the bandwidth of the 
search we can shorten the length of the data to be analyzed as we need to sample the 
data at only twice the bandwidth. For the sake of the EFT algorithm it is best to keep 
the length of the data to be a power of 2. Consequently we have chosen the number 
of data points to analyze to be A'' = 2-'^*. Thus for Tq = 2 days of observation time the 
bandwidth Ai/ was At/ — N / {2To) ^ 0.76 Hz. We have also chosen to analyze the data 
for the plus mode. As a result we searched the bandwidth from 921.00 Hz to 921.76 
Hz. We have used the filters with the phase linear in the parameters given by Eq. 
H17|l. In the filters we have included the amplitude modulation. The number of cells 
Nc calculated from Eq. 1)33(1 was around 1.6 x 10^^. Consequently from Eq. H34() the 
threshold signal-to-noise ratio for 1% false alarm probability was equal to 8.3. In the 
search that we have performed we have used a lower threshold signal-to-noise of 6.7. 
The aim of lowering the threshold was to make up for the loss of the signal-to-noise 
ratio due the discreteness of the grid of templates and due to the use of filters that 
only approximately matched the true signal. The number of points in the grid over 
which we had to calculate the statistics ^ turned out to be 183 064440. This number 
involved 63 830 positions in the sky and 2 868 spin down values for each sky position. 
We have carried out the search on a network of PCs and workstations. We had around 
two dozens of processors at our disposal. The data analysis started in September 2001 
and was completed in November 2002. 

The two-day stretch of data that we analyzed was taken from a larger set of 13 
days of data taken by the EXPLORER detector in November 1991. We have chosen 
the 2-day stretch of data on the basis of conformity of the data to the Gaussian random 
process. We have divided the data into 2^^ points sections corresponding to around 
11 hours of data. For each stretch we have performed the Kolmogorov-Smirnov (KS) 
test. The results of the KS test are presented in Figure Q For the case of KS test 
the higher the probability the more Gaussian the data are. From the above test we 

f The EXPLORER detector is operated by the ROG collaboration located in Italian Istituto 
Nazionale di Fisica Nucleare (INFN); seejhttp : //www ■ r omal . infn. it /rog/explorer/ explorer .html [. 
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Figure 1. Quality of the EXPLORER data. Tlie probability values of the KS 
statistics are given for each data segment from the 13 day data run. 



conclude that large parts of data are approximately Gaussian. On the basis of the 
above analysis we have chosen the two day data stretch to begin at the first sample of 
the 7th data stretch corresponding to Modified Julian Date of 48580.7909. In Figure|21 
we have presented the spectral density of the two day stretch of data that we analyzed. 
We see that the minimum spectral density was close to 10~^^/\/Hz. 

We have obtained 22 295 threshold crossings for the Northern hemisphere and 
44 701 for the Southern hemisphere. We consider all candidates contained within a 
single cell of the parameter space defined in Section 12.21 as dependent and we choose 
as an independent candidate the candidate that has the highest signal-to-noise ratio 
within one cell. Consequently we have obtained 11703 independent candidates for 
the Northern hemisphere and 18 702 for the Southern hemisphere. In Figures El and 01 
we have plotted the histograms of the values of the statistics J- for the independent 
candidates and we have compared it with the theoretical distribution for J- when no 
signal is present in the data. A good agreement with the theoretical distribution is 
another indication of Gaussianity of the data. It also reveals that there are no obvious 
populations of the continuous gravitational-wave sources at the level of sensitivity of 
our search. In the search no event has crossed our 99% confidence threshold of 8.3. 
The strongest signal obtained by the coarse search had the signal-to-noise ratio of 7.9 
and the fine search increased that value to 8.2. 

As we do not have a detection of a gravitational-wave signal we can make a 
statement about the upper bound for the gravitational-wave amplitude. To do this we 
take our strongest candidate of signal-to-noise ratio do and we suppose that it resulted 
from a gravitational- wave signal. Then, using formula (|24(l . we calculate the signal-to- 
noise dui of the gravitational-wave signal so that there is 1% probability that it crosses 
the threshold JFq corresponding to do, where JFq = 2-1- id^. The dui is the desired 99% 
confidence upper bound. For do — 8.2, which corresponds to the signal-to-noise ratio 
of our strongest candidate, we find that d^i = 5.9. For the EXPLORER detector this 
corresponds to the dimensionless amplitude of the gravitational-wave signal equal to 
2 X 10^^^. Thus we have the following conclusion: 

In the frequency band from 921.00 Hz to 921.76 Hz and for signals coming from 
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X 10" 



Srain sensitivity of the EXPLORER detector in the year 1991 




921.1 921.2 921.3 921.4 921.5 921.6 921.7 

Frequency [Hz] 



Figure 2. Two-sided spectral density of the EXPLORER data. 

Northern hemisphere 




Figure 3. Histogram of candidates: Northern hemisphere. On the x-axis the 
values of the statistic are given. The solid line corresponds to the theoretical 
distribution with 4 degrees of freedom. 



any sky direction the dimensionless amplitude of the gravitational-wave signal from a 
continuous source is less than 2 x 10~^^ with 99% confidence. 
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Southern hemisphere 
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Figure 4. Histogram of candidates: Southern hemisphere. 
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